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I'  The  problem  of  image  enhancement  by  nonlinear  two-dimensional  (2-D) 

1 

homomorphic  filtering  is  approached  using  stochastic  models  of  the  signal 
P and  degradations.  Homcmiorphlc  filtering  has  been  used  previously  for 

image  enhancement  but  the  lineeu*  filtering  operation  has  generally  been 

(i 

ii  chosen  heuristlcally.  In  this  paper  stochastic  image  models  described 

1 1 and  analyzed  previously  by  the  authors  are  used  to  model  the  true  image 

and  Interferrlng  components  (shadows  and  salt-and-pepper  noise).  The 
I problem  of  designing  the  linear  operation  ccua  then  be  formulated  as  one 

of  linear  least  mean-squared  error  (Wiener)  filtering.  Examples  of 
I processing  on  typical  real-world  images  are  Included  to  indicate  the 

, . results  obtainable. 
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I.  Introduction : 

A problem  in  representing  data  in  any  format  is  emphasizing  significant 
features  while  minimizing  distracting  detail.  Solutions  of  this  problem  when 
the  data  consists  of  images  are  called  image  enhancement  techniques. 

Elonents  in  Images  which  distract  the  viewer  from  relevant  Information 
are  shadows  and  salt-and-pepper  noise.  This  is  the  case  for  either  photo- 
chemical or  photoelectronlc  image  sensing  and  display  systems.  For  example, 
shadows  frequently  obscure  detail  in  ordinary  photographs;  even  when  deted.1 
is  present  within  the  area  darkened  by  the  shadow  it  is  often  subdued  because 
the  contract  in  the  darkened  curea  is  less  than  that  in  the  surrounding  bright 
regions.^ 

The  same  effect  occurs  in  the  recording  and  display  of  Images  by  photo- 
electronic  means  although  here  the  problem  is  further  aggravated  by  the  limit- 
ed dynamic  range  generally  available.  Many  TV  monitors  and  storage  CRT  dis- 
plays, for  exasqple,  aire  not  capable  of  isore  than  32  distinct  gray  levels.  The 
severity  of  this  restriction  is  beat  appreciated  in  comparison  to  the  capabil- 
ities afforded  by  a good  quality  photograph  under  subjective  humam  evaltiatlon. 
According  to  Zoethout  [l],  the  ratio  of  the  smeiLlesb  variation  in  intensity  dis- 
cernible by  the  eye  to  the  overall  intensity  is  generally  considered  to  be 
1/100.  Using  the  standard  definition  [2]  of  the  optical  density  D of  a 
photographic  emulsion 

D - logj^Q  1/T  , (1) 

where  T is  the  transmittance  of  the  emulsion,  we  conclude  that  the  smallest 
discernible  change  in  density  is  rou^ily  101/100b0.00U3.  laical  ^lotographic 

t This  is  due  to  the  common  practice  in  Bx>dem  photography  of  exposing  film 
so  that  shadows  fall  on  the  toe  of  the  D-log  E curve  (cf.  [19],  Fp.  29-32). 
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emulsions  (cf.  [2],  pg.  69)  cover  a range  of  densities  from  0 to  1.3> 

Hence,  assuming  grain  noise  Is  negligible  there  eure  approximately  3^7 
discernible  density  levels  In  a good  quality  black-and-vhlte  photograph; 
far  exceeding  the  capabilities  of  most  photoelectronlc  displays. 

The  mechcmlsm  by  which  the  eye  achieves  Its  vide  dynamic  range  suggests 
techniques  by  which  the  dynamic  range  of  photographic  materials  and  electronic 
displays  may  be  extended.  The  sensitivity  of  the  eye  to  local  changes  Is 
described  by  Hochberg  [3].  Within  a broad  range  of  illuminations  the  appeirent 
brightness  of  a region  Is  a function  of  the  ratio  between  Its  intensity  and 
that  of  Its  surrounding  region  and  not  of  its  absolute  intensity  alone.  It 
Is  as  if  the  eye  takes  the  logarithm  of  the  Incoming  light  intensities  and 
then  passes  these  signals  thro\i^  a spatial  highpass  filter  which  attenuates 
the  DC  ccaq>onent.  This  contention  Is  supi>orted  by  Stockham  [U]  who  carries 
It  to  the  point  of  estimating  the  frequency  response  of  the  highpass  filter. 
Additional  work  along  these  lines  has  teen  provided  by  Mannos  and  Sakrlson  [20] 

In  the  context  of  Image  coding. 

One  method,  then,  of  extending  the  dynamic  range  of  recorded  Images  is  to 
take  the  logarithm  of  the  original  intensities;  attenuate  slowly  varying 
components  which  contribute  little  information  to  the  eye  (and.  Indeed  act  as 
a source  of  degradation)  while  simultaneously  enhancing  rapidly  varying 
components  to  which  the  eye  is  most  sensitive;  and  finally,  presenting  the 
exponential  of  the  processed  image.  Photographers  have,  to  a degree,  acconq>llsh- 
ed  this  processing  in  recording  Isiages  on  film  by  using  a fast-acting  developer 
with  little  agitation;  exhausted  chemicals  from  heavily  exposed  areas  diffuse 
into  adjacent  lightly  exposed  areas  causing  them  to  develop  even  less  and 
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thereby  eiiq>ha8lzliig  transitions  from  light  to  daork  [2].  For  Images  stored 
in  digital  fonn  this  processing  is  readily  accomplished  by  a digital 
coBq>uter. 

This  paper  then  describes  a specific  digital  implementation  of  a scheme 
for  enhancing  rapid  local  variations  due  to  object  boundaries  and  de- 
ophasizlng  more  gradual  changes  such  as  those  which  occur  as  a result  of 
shadows  while  simultaneously  controlling  the  degrading  effects  due  to  salt- 
and-pepper  noise.  In  order  to  formulate  this  problem  the  observed  image 
intensity  ^^(x)  will  be  considered  to  be  the  product  of  three  caeq;>onents : 
the  true  reflectivity  of  the  objects  in  the  image  r(x),  the  non -constant 
illumination  function  l(x),  and  a white  noise  component  n(x).  That  is* 

Sq(x)  ■ r(x)  • i(x)  • n(x)  . (2) 

We  note  that  the  noise  here  is  modeled  as  a multiplicative  process  in  the 
intensity  domain  or,  equivalently,  additive  in  the  density  domain.  There  is 
substantial  evidence  to  Justify  this  for  either  photochemlceQ.  or  photo- 
electronic  sensing  and  display  systems.  The  dlscxission  in  Andrews  and  Hunt  [21] 
is  particularly  illunlnating  in  this  regard.  In  general,  however,  the  noise 
process  n(x)  will  depend  iqpon  the  signal  in  a rather  con^llcated  fashion.  We 
have  chosen  to  disregard  this  dependence  and  assume  a signal-independent  noise 
process  n(x).  Again  this  is  common  practice  in  dlgitcQ.  image  processing 
(cf.  [21],  pp.  20-23). 

Taking  the  logarithm  of  both  sides  of  (2)  then  res\ilts  in  the  density  image 
g(x)  - In  SQ(x)-f^(x)+fj^(x)+f^(x)  , (3) 

where  f (x)  ■ In  r(x),  f, (x)“  An  i(x),  and  finally  f (x)  ■ An  n(x).  If  precise 
r — n“  “ 

stochastic  descriptions  of  the  signals  f (x),  f. (x),  and  f (x)  are  available 

r “ 1 “ n “ 

the  problem  of  sorting  out  the  true  image  from  the  illumination  and  noise  can 
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be  formiilated  as  one  of  Wiener  or  linear  least  mean-square  filtering  where  the 
signed.  Is  f (x)  and  the  noise  Is  f. (x)+f  (x).  Passing  the  result  of  this  filter- 

I*  — 1 U 

Ing  through  an  exponentiation  operation  yields  an  estimate  f(x)  of  the  true 
Image  as  represented  by  the  reflectance  process  r(x).  A block  diagram  of 
the  Indicated  processing  Is  provided  In  Fig.  1. 

The  technique  of  sandwiching  a linear  operation  between  two  con^tllmentary 
nonlinear  operations  (homomorphic  filtering)  [7]  as  applied  to  Image  enhancement 
Is  not  new;  what  Is  \inlque  about  the  approach  described  here  Is  the  use  of 
explicit  stochastic  Image  models  as  a guide  In  the  choice  of  a linear  filter- 
ing oi>eratlon.  In  peortlcular,  we  develop  explicit  stochastic  models  for  both 
the  reflectance  process  f (x)  and  the  Illumination  process  f. (x).  The  reflect- 

ance  process  f (x)  hats  been  developed  to  exhibit  predominant  and  pronounced 
r — 

edge  structure  such  as  might  be  expected  of  real-world  Imagery.  This  process 

possesses  considerable  high  spatial  frequency  conqx^nents.  The  Illumination 

process  fj^(x),  on  the  other  hand,  has  been  specifically  developed  to  possess 

relatively  large  low  spatial  frequency  comi)onents  typical  of  Illumination  fields. 

These  stochastic  models  are  described  In  Sections  III  and  IV  respectively.  The 

noise  process  f (x)  Is  described  In  Section  V while  a description  of  the  digital 
n — 

Implementation  of  the  Wiener  filter  Is  provided  In  Section  VI.  Some  typical 
results  are  Illustrated  In  Section  VII  which  Is  followed  by  a summary  and  con- 
clusions In  Section  VIII.  We  begin  with  some  preliminary  discussions  on  random 
fields  as  Image  models. 

II.  Preliminary  Discussion : 

The  appro8M:h  described  eeurller  Is  not  complete  \intll  the  stochastic  de- 
scriptions of  f (x),  f. (z),  and  f (x)  are  specified.  The  models  to  be  considered 


p 

here  will  be  a family  of  random  variables  {f  (u),  xeR  } defined  on  the  plane 
where  u denotes  the  dependence  on  a fixed  probability  space  (fl,a,P).  Wong  [lO] 
refers  to  such  a collection  as  a random  field  but  notes  that  it  is  also  cailled 
a stochastic  process  with  a several-'dlmensional  time.  In  the  discussion  that 
follows,  for  convenience  the  dependence  upon  the  underlying  probability  space 
will  not  be  explicitly  denoted;  f(x)  will  be  written  for  Furthermore, 

consideration  will  be  restricted  to  zero-mean  fields  of  second  order  (veuri- 
ances  exist).  The  second-order  properties  are  described  by  the  covariance 
function  of  the  random  field  which  is  given  by 

= E {f(x)f(^)}  ; X,  ^eR^  , (U) 

where  E{ • > represents  the  expectation  operator  over  the  probability  space 
(n,a,P).  All  of  the  random  fields  to  be  presented  here  have  been  constructed 
so  as  to  possess  a covariance  function  invariant  \inder  all  Euclldecui  motions; 
according  to  the  definitions  of  Wong  [10 ] they  are  homogeneous  aind  isotropic. 
For  these  fields  (U)  can  be  written 

E{f(x+u)f(x)}=R^^(  I |u|  I );  u*^  = {\i^,  U2)eR2  , (5) 

where  | |u| | represents  the  Euclidean  norm  defined  in  terms  of  an  inner  product 
according  to  | |u| | * u^  + u|  . Since  the  random  fields  are  homogeneous 
their  power  spectreil  density  is  given  by 

Sff(i!»)  “ I Rff(u)exp  {-J<u,u>}du  , (6) 

r2 

T 

where  u ■ (u^,W2)  represents  a spatica  frequency  vector  emd  du  is  the 
differential  voliime  element  in  r2  . Fvurbhermore , for  isotropic  fields 
R^(u)»R^^(  I luj  I ) and  (6)  can  be  evaluated  up  to  functional  form  with  the 
aid  of  a theorem  of  Bochner  [ll]  yielding 


(7) 


Sj^(a))  » Sj.j.(n)  = 2irJ  XR^^(X)jQ(Xn)dX  , 

0 

where  (^1  |u|  I “ («i)|+o)|)  ^ represents  radial  frequency  and  Jq(*)  denotes  the 

ordinary  Bessel  function  of  the  first  kind  of  order  zero.  Thus  euid 

R^^(*)  are  related  throiigh  a Hankel  transform  [12],  [13]. 

The  enhancement  problem  as  ve  have  posed  it  consists  of  emphasizing 

f^(x)  while  minimizing  the  degrading  effects  of  f^Cx)  + specifically 

we  will  look  for  the  linear  least  meem-squared  estimate  of  f (x)  based  on 

r — 

observing  g(x).  If  we  have  knowledge  of  the  power  spectral  densities^ 

^d  associated  with  fj^(x),  and  ^^(x)  respectively, 

and  these  three  fields  are  uncorrelated,  the  2-dlmenslonal  (2-D)  linear  least 
mean-sq\uure  filter  possesses  transfer  function 


S ^(u) 

* Spr(i!L)S^i(]!>)+S^(«)  ’ 

which  is  simply  an  extension  to  2-dlfflenslons  of  the  Wiener  filter  [5]. 

We  will  now  proceed  tc  describe  the  random  fields  which  are  to  be  used  as 
the  reflectance.  Illuminance,  and  noise  models  respectively. 


III.  Reflectance  process ; 

The  process  to  be  considered  as  a model  of  the  reflectance  will  be  referred 
to  as  the  "rectangular  process".  This  process  divides  the  plane  into  regions 
whose  boundaries  form  rectwgles  €md  then  assigns  intensities  to  these  regions 
such  that  the  intensity  within  each  region  is  constcmt  and  possesses  specified 
correlation  properties  with  intensities  in  contiguous  regions.  It  was  originally 
developed  to  satisfy  the  need  for  an  edge  model  for  the  design  of  s linear  fllter- 

t In  general,  the  power  spectral  density  of  a random  field  f(x)  will  be  denoted 

by  ^ situations  where  subscripts  are  used,  such  as  f (x),  the 

eorrelpondlng  power  spectral  density  will  be  distinguished  by  the  associated 

subscript,  i.e.,  S^(u)  in  the  case  of  f (x). 

rr  “ r “ 
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ing  operation  to  be  used  in  edge  detection  Il5l.  Rougjily  speaking  it  is  a 
2-D  extension  of  a generalization  of  the  random  telegraph  vave  [lUJ, 


As  a first  step  in  describing  this  process, the  method  of  partioning  the 


pl€uie  into  rectangles  will  be  detailed.  A fundamental  role  in  this  stochastic 


image  model  will  be  played  by  the  Vector  valued  random  field  N(x)  which  prO' 


vides  a two-dimensional  generalization  of  a counting  process.  In  particular 


suppose  the  vector  x is  obtained  from  x according  to  x = ^ where  A is  the 


defined  for  some  0e[-Tr,v].  The  transformation  (9)  then  results  in  a rotation 
of  the  plane  in  terms  of  the  Cartesian  coordinate  axis  (x^ ,x„)  as  illustrated 


in  Fig.  2.  Consider  now  the  vector-valued  random  field  defined  by 


where  0 is  \miformly  distributed  on  [-it,ir]  and  {N^(t),  t ^ 0},  i=l,2  eure 
mutually  independent  one-dimensional  counting  processes,  i.e.,  N. (t)  repre^ 


sents  the  number  of  events  which  have  occurred  in  the  interval  [0,t]. 
Specifically, let  the  probability  density  function  p(t)  of  the  distance 
between  events)  be  given  by  the  exponentleil  distribution  function 


The  quantity  represents  the  average  number  of  events  per  \inlt  distance. 

Consider  now  the  random  field  (^x),  x^R^}  which  undergoes  transitions 
at  the  boundaries  of  the  elementary  rectangles  defined  by  {N(x),  xeR^}. 

The  gray  level  assumed  throughout  any  elementary  rectangle  will  be  assigned 


so  as  to  be  zei^s-mean  Gaussian'^with  vso’iance  and  correlated  with  the 

T 

gray  levels  in  contiguous  rectsuigles.  More  specifically,  assume  that  the 

random  orientation  0e [-»,»]  has  been  chosen  and  that  k transitions  have 

occurred  between  the  two  points  x and  x + u.  It  will  be  assumed  that 

E{f(3^u)f(x)l0,k>  = o^p^;  k=0,l,...,  (12) 

where  Ip  |<1  . The  quantity  p is  the  correlation  coefficient  governing  the 

spatial  evolution  of  the  random  an^lltude  process.  For  example,  let 

represent  the  amplitude  or  gray  level  assigned  after  i transitions  in  the 

direction  and  J transitions  in  the  x^  direction.  The  sequence{3^  then 

be  generated  recursively  according  to 

X = P X.  - ,+p  X.  . ,-P^X.  , , ,+W.  , , (13) 

i,j  r i-1,0  r 1,0-1  r i-i,o-i  i,0 

where  {W.  ,}  is  a 2-D  sequence  of  independent  and  identically  distributed 
i > J 

(l.i.d. ) zero  mean  Gaussian  variates  with  common  variance  o^=a^(l-p^)^. 

The  sequence  defined  by  (13)  has  an  alternative  interpretation  eis  the  output 

of  a recursive  2-D  digital  filter  excited  by  a white  noise  field. 

It  can  be  shown  that 

. 'Ik) 


so  that  the  condition  expressed  by  (12)  is  Indeed  satisfied. 

TypiceLl  computer-generated  realizations  of  the  resulting  random  field  are 
Illustrated  in  Fig.  3 for  selected  values  of  p^and  X^..  Note  that  X^ controls 
the  edge  density  while  p^  indicates  the  degree  of  correlation  of  intensities 
in  adjacent  regions.  The  displayed  Images  here  and  throughout  the  remainder 
of  this  paper  are  square  arrays  consisting  of  2^6  elements  or  samples  on  a 
side.  Vlhere  a value  of  X^is  specified  it  is  measured  in  normalized  units  of 

t The  Gaussian  assvmqrtlon  here  is  not  critical  and  is  easily  removed. 
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events  per  sample  distance  so  that  there  are  on  average  256  transitions  along 
each  of  the  orthogonal  axes. 

In  [15]  it  has  been  shown  that  the  power  spectral  density  of  this  process 
is  given  by 


^r(«)  = .2 


sd-ppxyoj 


2 I- 


n"+2(l-p^2x2 


(15) 


Although  realizations  of  this  process  may  not  resemble  any  particular  real- 
world  image, the  general  properties  of  regions  of  constant  reflectivity  and 
boundaries  of  regions  being  parallel  are  common  to  many  real-world  images. 

More  specifically,  we  do  not  claim  that  the  model  proposed  here  is  universally 
representative  of  real-world  imagery  except  in  a certain  qualitative  sense. 

This  model  is  completely  defined,  up  to  a scale  factor,  by  the 
two  peirameters  X^  and  p^.  The  parameter  X^  represents  the  "edge  busyness" 
associated  with  an  image  while  p^  is  indicative  of  the  degree  of  abruptness 
across  an  edge  boundary.  For  p^  large  (in  magnitude)  and  negative  there  is 
an  abrupt  almost  black-to— white  or  white-to-black  tremsition  across  an  edge 


boundary.  For  p>0,  on  the  other  hand,  the  transition  across  an  edge  boimdary 
is  much  more  gradual.  These  properties  are  easily  related  to  any  particular 
real-world  image. 

IV.  Illumination  Process; 

In  the  previoTiB  section  a 2-D  random  field  which  possesses  em  inherent 
rectangular  structvire  was  described  and  analyzed.  While  this  random  field 
is  a reasonable  model  of  objects,  shadows  do  not  generally  exhibit  this  rec- 
tangular mosaic  but  rather  a much  more  random  edge  orientation.  Here  the  con- 
struction and  properties  of  a class  of  2-D  random  fields  which  provide  a more 
appropriate  model  for  lUumlnation  in  the  sense  that  individual  realizations 
do  not  exhibit  this  rectangular  pattern  will  be  presented.  Consider  the  random 
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partition  of  the  plane  by  a series  of  sensed  (directed)  lines  which  can  be 
described  as  a marked  point  process.  Specifically,  an  arbltraiy  sensed  line 
can  be  described  in  terms  of  a 3-tuple  (r,9,d)  where  r represents  the 
perpendicular  distance  or  radius  to  the  line  in  question,  6e[-ir  ,ii  ] 
represents  the  orientation  of  the  vector  of  lent^th  r drawn  from  the  origin 
and  perpendicular  to  this  line,  and  de{-l,  1}  specifies  the  sense  of  the 
line  where  Fig.  U illustrates  a case  of  d=l  (pointing  counterclockwise). 

As  shown  in  Fig.  U,  by  virtue  of  the  direction  imposed  on  this  line  segment 
the  plane  is  partitioned  into  two  distinct  regions,  R(rlght  of  line)  and 
L(left  of  line)  such  that  RUL*R^* 

Now  consider  the  set  of  lines  represented  by  the  sequence  {(rj^,0j ,d^)}. 
Here  {r^^}  represents  the  event  positions  associated  with  a homogeneous 
Poisson  process  {N(r),r  ^ 0}  with  intensity  X events/vmit  distance  evolving 
according  to  the  radial  parameter  r.  The  sequence  {P^}  represents  an 
(i.i.d. ) sequence  of  random  variables  uniformly  distributed  on  [-Tr,ir] 
while  the  sequence  {d^^}  is  likewise  i.i.d.  assuming  the  binary  values  ±1 
with  equal  probability.  Finally,  the  sequences  {r^^},  {6^}  and  {d^^} 
are  assumed  to  be  independent. 

Now  note  that  the  set  of  lines  generated  by  the  process  described  above 
will  divide  the  plane  into  disjoint  polygonal  regions.  For  any  point  x 
define  Njj^(3c)  to  be  the  number  of  left-to-right  crossings  of  sensed  liner 
incurred  in  moving  along  a stredght  path  from  the  origin  to  the  point  x; 
similarly  define  Hpj^(x)  to  be  the  number  of  right-to-left  crossings.  Observe 
that  Kj^(ij^)  * ^LR^ig^  ®RL^-1^  * points  x^  and  Xg 

which  are  in  the  same  region.  Assign  a label  k to  each  region  as  k » Nj^(x)  - 
Nj^(x)  where  x is  a point  in  the  region.  This  label  can  in  turn  be  used  to 
assign  a gray  level  to  each  region  by  using  it  to  index  into  a sequence  (X^^) 
of  zero-mean  Gaussian  random  variables.  In  particular,  assume  that  the 
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sequence  {x^^}  is  generated  recursively  according  to 

\ = piVi  ^ ''k  • 


»k*l,2,...  f 


\ ’ PiVl  " \ 


} k * —1  ,—2 1 , 


vitb  {Wj^}  cm  l.i.d.  zero  mean  Gaussian  sequence  vlth  variance  « (l-Ppo^ 
and  the  initial  value  chosen  to  have  zero-mean  Gaussian  distribution  vith 
variance  o^.  Note  that  this  sequence  is  stationary  with  covcuriance  given  by 

Selected  computer-generated  realizations  of  the  resulting  random  field 
are  Illustrated  in  Fig.  5 for  various  values  of  emd  * X/ir.  The  quantity 
X^  can  be  shovn  to  represent  the  average  edge  density  along  any  randomly 
chosen  line  segment  on  the  plane.  From  Fig.  3 observe  that  for  p - 0.3  light 
regions  tend  to  be  surrounded  by  light  regions  while  for  p^^  = -0.9  there  is 
a tendency  for  light  regions  to  be  surrounded  by  dark  regions  and  vice  versa. 
Indeed,  in  the  latter  case  there  is  almost  a black-to-white  or  white-to-black 
Inversion  across  an  edge  bo\indary. 

An  expression  for  the  autocorrelation  function  given  by  (U)  ccm  be  obtained 
as  in  [17]  with  the  result 

Rii(u)  * a*e  {Iq(Xj|1u||)  + 2 ^ |u|  | ) ) , (18) 

k— 1 

where  is  the  modified  Bessel  function  of  the  first  kind  of  order  k. 

Simllcurly,  the  corresponding  power  spectrcuL  density  evaluated  according  to  (7) 
is  given  by 


Su(n). 


f»r  1. 


■*’i'  P r 1-C08»  dj ^ 

I l-2p^C08t+pJ  [(n/>  )2+(i_eo8*)2]^2 

0 L 
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One  noteable 


which  is  Illustrated  in  Fig.  6 for  various  values  of 

characteristic  of  this  random  field  is  that  for  Q small  compared  to 

—3/2 

the  power  spectral  density  behaves  approximately  as  fl  ' , i.e., 
has  a singularity  at  the  origin  except  for  = -1.  This  high  concentration 
of  energy  at  low  spatial  frequencies  is  a direct  result  of  the  construction 
procedure  which  allows  relatively  large  correlation  between  gray  levels  in 
regions  relatively  far  apart . We  feel  that  this  characteristic  is  typical 
of  random  illumination  processes  and  as  a result  it  was  purposely  built  into 
the  construction  procedure. 

The  effect  on  a real- woritt  image  of  multiplying  the  intensity  by  the  ex- 
ponentieil  of  a realization  of  the  polygonal  process  is  Illustrated  in  Fig.  7. 

Not^  how  the  illumination  process  seems  simply  to  have  added  extra  shadows  to 
the  building.  Although  somewhat  of  an  exageration,  this  figure  has  been  includ- 
ed to  illustrate  the  nature  of  effects  which  may  be  caused  by  the  proposed  random 
illumination  process. 

V.  White  Noise  Process; 

In  any  image  that  is  scanned  there  is  some  noise  present  that  is  independ- 
ent from  pixel-to-pixel;  due  to  its  appeaurance  it  is  often  referred  to  as 
"salt-and-pepper  noise”.  Here  we  will  model  this  noise,  which  is  represented 
by  as  white  Gaussian  noise.  While  the  Gaussieui  assumption  is  not  critlcail 

to  the  analysis  which  follows  it  is  interesting  to  note  the  observation  of  Selwyn 
(cf.  [2],  pg.  12^)  on  noise  due  to  film  grains.  He  notes  that  if  the  spot  area 
covered  by  a beam  used  to  measure  the  density  of  a photographic  em\ilsion  is 
much  larger  than  the  projection  of  individual  gredns  then  the  density  distri- 
bution will  be  Gaussian  if  the  grain  noise  is  small,  and  presumably  homogeneous. 
In  fact,  the  following  relationship  between  the  size  of  the  scanning  aperture 
and  variance  of  the  measured  densities  has  been  found  to  hold: 


o2  = G2/2a  , 
n 


(20) 


where  is  the  variance  of  the  density,  G is  a measure  of  the  granularity 
of  the  film,  and  a is  the  area  of  the  scanning  aperture.  If  the  scanning 
aperture  is  considered  to  be  a spatial  lowpass  filter  whose  bandwidth  is 
proportional  to  l/v^  then  (20)  implies  that  the  noise  power  present  in  the 
density  measurement  is  proportional  to  the  lowpass  filter  bandwidth  squared. 
Since  this  is  the  effect  that  would  be  observed  if  the  grain  noise  were  white 
this  modeling  assumption  is  Justified.  The  power  spectral  density  of  this  noise 
can  be  expressed  then  as 

S„„(n)  = a2  (21) 

nn  n ’ 

where  is  the  variance  of  the  samples. 

VI.  Linear  Filter  Description  and  Implementation; 


Using  the  models  presented  above  the  linear  least  mean-squared  error 
(Wiener)  filter  for  estimating  the  reflectemce  process  will  now  be  derived. 
The  system  transfer  function  given  by  (8)  cam  be  rewritten  as 


S ( w ) S_^(w) 


(22) 


rr— • rr'— ’ 

which  suggests  defining  the  parameters  and  as  the  ratics  of 
illumination  and  noise  powers  to  reflectance  power  respectively.  These 
parameters  provide  a measure  of  the  degree  of  degradation  caused  by  shadows 
and  salt -and -pepper  noise.  More  specifically,  we  have 


while 


i 

(23) 

(2U) 
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so  that  the  filter  is  oonpletely  specified  In  terms  of  the  quantities 

Yi»  Y » Py*  Finally  .observe  that  for  the  given  power  spectral 

densities  the  filter  is  radially  symmetric,  i.e.,  it  is  a function  ofn»|lu||. 

Since  the  final  objective  of  this  process  is  to  make  the  best  use  of  the 
limited  dynamic  range  avcdlable  for  display  it  will  prove  expedient  to  include 
a gain  factor  chosen  such  that  the  power  leaving  the  linear  processor  equals 
that  entering.  it  follows  then  that  is  given  by 

Aq  = ///  dw  , (25) 


where  S (u)  is  the  power  spectral  density  of  the  filter  input  gCx)  in  (2)  ahd 


Is  given  by 


gg 


li' 


nn 


(26) 


The  description  of  the  filter  provided  above  is  for  a continuous  spatial 
domain;  ve  have  access  only  to  ssunpled  data.  Hence  we  seek  a C2-D) 
digitfd.  filter  with  system  transfer  function  whose  frequency  response^ 

approximates  HQ(n).  Specifically , the  Wiener  filter  has  been  implemented  as  a 
2-D  infinite  impulse  response  (HR)  [18]  digital  filter  whose  point 
spread  function  exhibits  quadrant  symmetry.  The  choice  of  an  HR  filter  was 
based  on  the  computational  economies  which  result  from  the  ability  to  imple- 
ment the  filter  with  a recursive  structure.  Consider  first  a filter  with  the 
following  recursive  structvire: 


* ’ (1  j)eD' 


The  frequency  response  of  the  tyo-dimensional  digital  filter  is  simply 
®0^*1**2^  evaluated  for  z^*  eJ“i,  1»1,2.  We  assume  the  spatial  frequency 

variable  is  measured  in  units  of  radians  per  sample  distance. 

ll( 


where^  } and  D'  * D - {(0,0)}.  At  this 

point  ve  will  assme  that  the  coefficients  of  y in  (32)  have  been  chosen  such 
that  the  resulting  filter  is  stable  when  recursing  from  the  iq>per  left-hand 
comer.  Observe  that  the  resulting  filter  vlll  have  nonzero  response  only  in  the 
lower  right  quadrant.  Since  the  filter  described  in  (22)  exhibits  radial  symmetry 
it  is  desirable  for  the  resulting  digital  filter  to  exhibit  four-quadrant  symmetry. 

One  method  for  achieving  a point  spread  function  with  this  Inherent 
symmetry  is  to  allow  repeated  application  of  the  same  filter  recursing  from 
each  of  the  four  corners.  If  {h.  ,}  represents  the  point  spread  function 
associated  with  a single  application  of  the  filter  8i>eclfied  by  (27)  then  the 
composite  filter  possesses  point  spread  function  as  indicated  in  Fig.  8.  The 
corresponding  system  transfer  function  of  the  composite  filter  is  then 

= Hq(z^.Z2)+Hq(z]^^  Z2)+Hq(z^,z’^)+Hq(z"^,z“^),  (28) 

where  is  the  transfer  function  corresponding  to  (27)  and  is  given  by 


° ^ 2 (i,j)eD  ^ 


In  choosing  li^{z^,z^)  to  provide  an  approximation  to  8^(0)  for  z^=e'^‘*^  , i=l,2 
we  have  restricted  attention  to  the  case  where  Hq(Zj^,Z2)  is  a second  order 


section  of  the  form 


Ho(Zi.Z2) 


^%1^  ^^^11  ^ *1^*2^  ^%2<  ^ 


with  the  restriction 


^00*  “ ^01 


-\l- 


This  choice  Insures  zero  frequency  response  at  the  origin  and  symmetry  of  the 


The  8iq>port  D of  the  recursive  2-D  filters  taken  here  will  be  such  as  to 
force  a quartet*  plane  causality  condition.  That  is,  D consists  of  a finite 
number  of  points  above  and  to  the  left  of  the  pixel  in  position  (i.j). 
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corresponding  point  spread  function  about  a line  U5®  to  the  axis;  l.e.,  h.  =h  .. 

i > J J >* 

Similar  properties  extend,  of  course,  to  the  composite  filter  represented  by 
H(z^,Z2).  a computer  program  has  been  written  for  determining  the  coefficients 

'*00’  ^01’  ^11*  ^02’  ®'01’  *^11  *^02  'to  an  Iterative  gradient  proced- 

ure to  result  In  a frequency  response  for  H{z^,Z2)  which  provides  a least  mean- 
squared  error  approximation  to  the  desired  response  11^(0).  The  details  of  this 
program  are  described  In  [6].  In  Table  1 we  summarize  the  results  of  this 
procedure  for  selected  values  of  P^»  and  p^,  all  with  Yp=-20dB. 

Here  we  have  foxind  It  convenient  to  classify  the  shadow  characteristic  as  soft, 
medium  or  harsh  depending  upon  the  value  of  Yj^* 


Shadov 

Characteristic 


Soft 

Medlus 

Harsh 


.0125 

0.5 

.05 

.0125 

0.5 

.05 

.0125 

-.9 

.05 

I 


*01 

“ll 

*02 

'*00 

\l 

^1 

'*02 

-.8650 

.7830 

-.oaUo 

.2376 

-.1957 

.1831 

-.0lU7 

-.8503 

.77U3 

-.0326 

.2661 

-.2216 

.2072 

-.0160 

-.8817 

.7959 

-.0118 

.5109 

-.U718 

.Ult58 

-.0065 

Table  1 

lypical  Filter  Parameters  for  Yj^*“20dB 


For  the  three  cases  described  in  Table  1,  the  corresponding  power  spectral 
densities  S^(n)  imd  together  with  the  resulting  Wiener  filter  response 

HqCJI)  are  plotted  in  Fig.  9 as  a function  of  the  radial  frequency  variable  il. 
Observe  the  sharp  DC  rejection  in  all  cases.  This  Is,  of  course,  due  to  the 
relatively  high  concentration  of  energy  at  low  spatial  frequencies  associated 
with  the  illumination  process  t^ix)  as  described  previously.  The  Wiener  filter 
response  Hq(0)  exhibits  the  essential  characteristics  of  the  linear  hlghpass 
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filter  originally  proposed  by  Stockham  on  an  ad  hoc  basis.  It  is  of 
some  interest  that  this  choice  can  be  Justified  rigorously  under  the  model- 
ing assumptions  made  here. 


In  Fig.  10  we  Illustrate  3-dimensional  plots  of  one  quarter-plane  of  the 
frequency  response  of  the  filters  associated  with  the  three  cases  described 
in  Table  1.  The  left  hand  column  shows  the  desired  responses  while  the  right 
hand  column  illustrates  the  responses  obtained  by  the  digital  filters.  The 
frequencies  displayed  run  from  0 at  the  front  to  ir  radians  per  pixel  at  the 
back.  Note  that  in  all  cases  the  complete  atten\iation  of  DC  and  rolloff  at 
higher  frequencies  have  been  preserved  in  the  transition  from  the  desired 
analog  filters  to  the  digital  approximations. 


VII.  Results: 

The  algorithm  presented  above  has  been  applied  to  several  real  world  images 
for  selected  values  of  the  parameters  p^,  X^,  p^,  and  X^.  The  processed 

Images  shown  here  Illustrate  the  effects  of  varying  these  parameters  on  the 
reconstructed  Images. 

Parts  a of  each  of  the  foiir  sets  of  images  composing  Fig.’s  11  through  lU 
are  the  original  Images  while  parts  b,  c,  and  d are  the  corresponding  processed 
versions.  All  of  the  processed  images  have  the  salt-and-pepper  noise  20  dB 
below  the  reflectance  and  the  reflectance  has  zero  correlation  of  intensities 
in  different  regions  (p  *0)  and  em  edge  density  such  that  on  the  avereige  there 

T 

will  be  12  transitions  on  a side  of  the  image  ( X ^.05).  The  variations  among 

r 

b,  c,  and  d are  due  to  chemges  in  the  illumination  intensity  y^^*  and  the 
illumination  correlation  p^. 
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Part  b of  these  figiires  has  the  edge  density  of  the  illumination  process 
= .0125  or, on  the  average, 3 transitions  on  a side  of  the  image,  p^»0.5 
(positive  correlation  of  illumination  in  adjacent  regions)  and  illumination 
power  6 dB  below  reflectance  power  dB),  corresponding  to  a soft  shadow. 

Note  how  in  part  b local  variations  are  much  more  evident;  shadow  and 
highlight  detail  have  both  been  emphasized.  Specifically,  in  Fig.  Ub  detail 
in  the  hair  is  much  more  evident  than  in  lla.  In  Fig.  12b  the  light  scaffolding 
blends  in  less  with  the  clouds  while  the  dark  slatting  on  the  side  of  the  tower 
is  also  more  evident.  Similarly,  in  Fig.  I3b  the  bri^t  regions  in  the  center 
have  been  de— emphasized  and  detail  in  the  building  is  more  apparent.  Finally , 
for  the  forward  looking  infrared  (FLIR)  Image  in  Fig.  lUb  the  tank  blends  in 
less  with  its  soirroundings . This  processing  brings  out  detail  in  both  dark 
and  light  areas. 

Parts  c (medium  shadow)  of  Fig.'s  11-lU  vary  from  parts  'b  in  the  increase 
of  the  illumination  power  by  3dB.  The  effects  noted  above  have  simply  become 
more  evident. 

Part  d (harsh  shadow)  deviates  from  c in  both  illumination  power  (which  is 
increased  by  9dB)and  the  correlation  of  illuminances  in  different  regions  which 
is  now  -0.9  instead  of  0.5.  Part  d of  all  the  figures  brings  out  detail  in 
adjacent  light  and  dark  regions  better  them  c as  is  expected  from  the  shadow 
being  modeled  as  having  balck-to-white  transitions  cu:ross  boundaries  between 
adjacent  regions. 


j 
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S 
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An  application  of  this  technique  to  a dynamic  range  reduction  problem  Is 
Illustrated  In  Fig.  1^.  In  the  original  tomograph  Image,  little  detail  Is 
evident  In  the  dark  area  while  the  light  areas  appear  muddy.  In  Fig.  l^b  some 
detail  Is  evident  In  the  dark  region  while  local  shading  Is  more  appaurent  In 
light  regions.  VHien  harsher  shadows  are  assumed,  as  In  l^c  emd  l^d,  local 
detail  Is  amphaslzed. 

In  all  cases  the  processing  has  Increased  the  visibility  of  detail  In  such 
a way  that  less  scrutiny  need  be  applied  when  observing  the  processed  imeiges 
than  is  required  when  observing  the  originals.  There  are,  however,  differ- 
ences in  the  chauracteristics  of  the  processing  that  occurs.  Increasing  the 
shadow  power  auid  making  the  shadow  more  pronotuiced  by  having  negative  correl- 
ation (p*-0.9)  instead  of  positive  correlation  (p=0.5)  has  made  local  detail 
more  visible  by  suppressing  variations  in  Intensity  that  occur  over  large  regions. 
VIII  Simmary  and  Conclusions: 

A method  for  determining  the  linear  filtering  operation  to  be  performed 
as  paurt  of  homomorphic  processing  based  on  Wiener  filtering  concepts  has  been 
described.  The  resulting  linear  filter  is  determined  by  the  parameters  of  the 
assumed  stochastic  models  of  the  signal  auid  degradation.  These  parameters 
aure  associated  with  salient  characteristics  of  the  signal  auid  degradation. 
Although  no  claims  cam  be  made  about  the  optimaaity  of  a particular  filter 
for  any  given  real-world  Imaige  It  has  been  shown  that  the  parameters  in  the 
models  provide  a useful  meams  for  specifying  the  desired  effects  of  the 
processing. 

The  Implementation  of  the  linear  operation  has  tadeen  the  form  of  a 2-D  re- 
cursive IIR  digital  filter  whose  order  ham  been  restricted  to  result  in  compu- 
tationad  efficiency  amd  whose  coefficients  have  been  chosen  to  result  in  a good 
approximation  to  the  Wiener  filter. 
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The  resulting  nonlinear  filtering  operation  has  been  applied  to  a 
number  of  real-vorld  images  to  illustrate  the  relationship  between  the 
parameters  of  the  assumed  stochastic  models  and  the  effects  of  the  re- 


References 


1.  W.D.t  Zoethout  , Physiological  Optics.  The  Professional  Press,  Chicago,  111. 
pp.  157,  1939. 

2.  P.  Kowaliski  , Applied  Photographic  Theory.  Wiley,  New  York,  1972. 

3.  J.  E.  Hochberg  , Perception.  Prentice-Hall,  Englewood  Cliffs,  N. J. , I96U. 

U.  T.  G.  StOCkham  , "Image  Processing  in  the  Context  of  a Visual  Model", 

Proceedings  of  the  IEEE,  Vol.  60,  No.  7,  pp  828-8U2,  July  1972. 

5.  H.  L.  Van  Trees  , Detection.  Estimation  emd  Modulation  Theory;Vol.  1.  Wiley, N.Y. 1968. 

6.  R.  W.  Fries  , "Theory  and  Application  of  a Class  of  Two-Dimensional  Random 

Fields",  Ph.D.  Thesis,  RPI,  Troy,  N.Y. , in  preparation. 

7.  A.  V.  Oppenheim,  R.  W.  Schafer  ,and  T.G.  Stockham,  Jr.,  "Nonlinear  Filtering  of 

Multiplied  and  Convolved  Signals",  Proc.  IEEE,  Vol.  56,  pp  126U-1291,  Aug.  I968. 

8.  E.  Wong  , "Homogeneous  Gauss-Maurkov  Rauidom  Fields",  Ann.  Math.  Stat.,Vol.  Uo, 

pp.  I625-I63U,  1969. 

9«  E.  Wong  , "Two-Dimensional  Random  Fields  and  the  Representation  of  Images", 

SIAM  J.  Appl.  Math., Vol.  16,  pp.  756-770,  1968. 


1 

li 

0 

n 


10.  E.  Wong  , Stochastic  Processes  in  Information  and  Dynamical  Systems.  Chap  7, 

McGraw-Hill,  New  York,  1971. 

11.  S.  Bochner  , Lectures  on  Fourier  Integrals.  Annals,  of  Math.  Studies,  Noi  h2, 

Princeton  Univ.  Press,  Princeton,  N.J.,  pp. 235-238,  1959. 

12.  A.  Papoulis  , "Optical.  Systems,  Singiilarity  Functions,  Complex  Hankel  Transforms", 

J.  Opt.  Soc.  Amer.,  Vol.  57,  pp*  207-213,  1967. 

13.  A.  Papoulis  , Systems  auad  Transforms  with  Applications  in  Optics.  McGraw- 

Hill,  New  York,  1968. 

1^*.  J.  W.  Modestino  auid  R.  W.  Fries  , "A  Generailization  of  the  Random  Telegraph 
Wave",  submitted  to  IEEE  Trauis.  on  Inform.  Theory. 

15.  J.  W.  Modestino  and  R.  W.  Fries  , "Edge  Detection  in  Noisy  Images  Using 

Recursive  Digital  Filtering",  Computer  Graphics  amd  Image  Processing,^, 
pp.  U09-U33,  Oct.  1977. 

16.  J.  W.  Modestino  , R.  W.  Fries  and  0.  G.  Daut  , "A  Generaaization  of  the 

Two-Dlmensionad.  Rimdom  Checkerboaurd  Process", to  appeau:  in  J.  Opt.  Soc.  Amer. 

17.  J.  W.  Modestino  amd  R.  W.  Fries  , "Construction  amd  Properties  of  a Useful 

Two-Dimenslonail  Random  Field",  submitted  to  IEEE  Trams,  on  Inform. Theory. 

18.  L.  R.  Rablner  and  B.  Gold,  Theory  and  Applications  of  Digital  Signal  Processing. 

Prentice-Hall,  Englewood  Cliffs,  H.J.,Pg.U51,  1975* 


21 


k. 


19.  C.  Shipman,  Understanding  Photograi 


H.  P.  Books,  Tucson,  Arizona,  197** 


20.  J.  Mannos  and  D.  L.  Sakr Ison, "The  Effects  of  a Visual  Fidelity  Criterion 
on  the  Encoding  of  Images",  IEEE  Trans.  Inform.  Theory,  Vol.  20,  pp.  525' 
536,  197**. 


21.  H.  C.  Andrews  and  B.  R.  Hunt,  Digital 
Englewood  Cliffs,  HJ,  1977- 


:e  Restoration.  Prentice-Hall 


Figure  3 

Selected  Realizations  of  Reflectance  Process 
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a.)  Original 


b.)  Illumination 


Original  with  Illumination 


Figure  7 


Effect  of  Illumination  Process 
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Figure  8 

Point  Spread  Function  of  Composite  Filter 
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a.)  Desired  Response  for  b.)  Actual  Response  for 
Soft  Shadow  Soft  Shadow 


c.)  Desired  Response  for  d.)  Actual  Response  for 
Medium  Shadow  Medium  Shadow 


e.)  Desired  Response  for  f.)  Actual  Response  for 
Harsh  Shadow  Harsh  Shadow 


Figure  10 

Frequency  Responses  of  Desired  and  Actual  Filters 


c . ) Medium  Shadow 


d. ) Harsh  Shadow 


Figure  11 

Typical  Results  on  Head-and-Shoulder  Image 
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a. ) Original 
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Figure  13 
Typical  Results  on  Building  Image 
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b . ) Soft  Shadow 


c.)  Medium  Shadow 


Typical  Results  on  Tomograph 
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